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Abstract 

Delamination growth caused by local buckling 
of a delaminated group of plies was investigated. 
Delamination growth was assumed to be governed 
by the strain-energy release rates G/, G//, and G///. 
The strain-energy release rates were calculated using 
a geometrically nonlinear, three-dimensional, finite 
element analysis. The program is described and 
several checks of the analysis are discussed. Based on 
a limited parametric study, the following conclusions 
were reached: 

1. The problem is definitely mixed-mode. In 
some cases G/ is larger than G//; for other 
cases the opposite is true. 

2. In general, there is a large gradient in the 
strain-energy release rates along the delami- 
nation front. 

3. The locations of maximum G/ and Gjj de- 
pend on the delamination shape and the 
applied strain. 

4. The mode III component was negligible for all 
cases considered. 

5. The analysis predicted that parts of the delam- 
ination would overlap. The results presented 
herein did not impose contact constraints to 
prevent overlapping. Further work is needed 
to determine the effects of allowing the over- 
lapping. 

Introduction 

Composite materials exhibit a number of different 
failure mechanisms. Among these are delamination, 
fiber breakage, and intralaminar cracking. In addi- 
tion to the failure mechanisms operative under ten- 
sile loads, there are mechanisms related to instability, 
such as fiber microbuckling and buckling of a group 
of delaminated laminae, which can be important un- 
der compression loads. Delamination which is driven 
by local buckling of a group of laminae is discussed 
in this report. This mechanism is referred to herein 
as instability-related delamination growth (IRDG). 
The magnitudes of the modes of strain-energy release 
rates G/, Gji, and G/// are often used to predict 
when delamination will occur. It is assumed that 
G/, G//, and G/// are the critical parameters for 
predicting delamination growth. 

Figure 1 shows two configurations which exhibit 
IRDG, laminates with a through-width or an embed- 
ded delamination. The through-width delamination 
has been analyzed using a variety of detailed and ap- 
proximate stress analyses (refs. 1 to 6). Experimen- 
tal measurements of IRDG have been published for 


static loads (refs. 6 to 9) and fatigue loads (refs. 1, 
3, 7, and 8). The through-width delamination is per- 
haps the simplest configuration that exhibits IRDG. 
Hence, it provides a convenient vehicle for evaluating 
various ideas about testing and analysis. 

The embedded-delamination configuration has 
also been studied both analytically (refs. 10 to 15) 
and experimentally (refs. 7, 12, and 16 to 18). To 
date, however, there have been no detailed analyses 
of the embedded delamination. The analyses have 
been limited to Rayleigh-Ritz or finite element plate 
analysis. With one exception (ref. 14), only aver- 
age total strain-energy release rates along the delam- 
ination front have been presented. In reference 14, 
the distribution of total strain-energy release rate 
was calculated. A geometrically nonlinear, three- 
dimensional, finite element analysis (or some other 
numerical technique) is required to calculate the in- 
dividual modes of strain-energy release rate. How- 
ever, three-dimensional finite element analysis tends 
to be very expensive. Since nonlinear analysis re- 
quires iteration, it is even more expensive. This ex- 
pense probably explains the lack of detailed analysis 
for this configuration in the literature. However, un- 
til detailed analyses are performed, the accuracy of 
approximate analyses cannot be assessed. 

There are two objectives to this paper. The first is 
to describe the three-dimensional, geometrically non- 
linear, finite element program which was developed 
for this study. The program is named NONLIN3D. 
The second objective is to show the effect of several 
parameters on the magnitude and distribution of the 
strain-energy release rates. The parameters consid- 
ered were delamination shape, delamination size, and 
strain level. 

Symbols 

a, b 


Cii dii 




^ 22 ) #33 

pa 

G 


semi-axes of elliptical 
delamination in x- and 
^-directions, respectively 

nodes used in calculation 
of strain-energy release 
rates 

constitutive coefficients 

constitutive coefficients 
for a “homogeneous 
quasi-isotropic laminate” 

Young’s moduli for 
orthotropic material 

nodal forces 

strain-energy release rate 



Gi, Gji, 

G///» Gt 

mode I, mode II, mode 
III, and total strain- 
energy release rates 

G 12, G 23 

, G 13 

shear moduli for or- 
thotropic material 

H 


thickness of base 
laminate 

h 


thickness of sublaminate 

K *P 


coefficients in tangential 
stiffness matrix 

P 


transverse load 

q a 


nodal displacements 

S 


perimeter coordinate 

u 


strain energy 

u, v, w 


displacements in x-, y-, 
and ^-directions 

V 


volume 

W 


width of finite element 
model 

W 0 


transverse displacement 
in center of plate 

X, y , 2 


rectangular Cartesian 
coordinates in global 
coordinate system 

x', y', z' 


rectangular Cartesian 
coordinates in local 
coordinate system 

A a 


increment in delamina- 
tion length 



strains 

£x 


strain in ^-direction 

V 12 , ^23) 

v 13 

Poisson’s ratios for 
orthotropic material 

n 


total potential energy 


Analysis 


The theoretical aspects of NONLIN3D, the ge- 
ometrically nonlinear, three-dimensional, finite ele- 
ment program used for this study, are discussed in 
this section. In particular, the following topics are 
covered: 


1. Governing nonlinear equations 

2. Substructuring 

3. Calculation of strain-energy release rate 

After the discussion of NONLIN3D, the mesh gen- 
eration and the finite element models are discussed. 
The material properties are also described. 

Governing Nonlinear Equations 

The derivation of the equilibrium equations and 
the expressions for the internally generated nodal 
forces and the tangent stiffness matrix are discussed 
in this subsection. The strain-displacement relations 
are also discussed. 

The total potential energy is given by 

U = \J CijSiSjdV - F a q a (1) 

where the integral term is the strain energy and the 
second term is the potential energy of the applied 
loads. The terms CC- and £j are terms in the constitu- 
tive matrix and the strains, respectively. The terms 
F a and q a are the generalized forces and displace- 
ments, respectively. The adjective “generalized” is 
used to indicate that F a and q a need not be forces 
and displacements in the usual sense. For example, 
in traditional Rayleigh-Ritz analyses, the unknowns 
q a are simply coefficients in the series expansion for 
the displacements. In this discussion, however, F a 
and q a always refer to the nodal forces and displace- 
ments in the x-, y-, and ^-directions. The system is 
assumed to be conservative; hence, the equilibrium 
state is obtained by minimizing n, which is accom- 
plished by setting the first partial derivatives with 
respect to the unknowns equal to zero: 

|2- = f CijSi dV -F a = 0 (2) 

d q a J V l dq a v ' 

Equation (2) is nonlinear because of the nonlinear 
strain-displacement relations. The integral in equa- 
tion (2) gives the magnitude of the internally gener- 
ated nodal forces corresponding to the current dis- 
placements. Until a converged solution is obtained, 
the internally generated forces do not equal the ex- 
ternally applied forces. The differences in the forces, 
referred to as residuals, equal The Newton- 

Raphson procedure that was used to solve equa- 
tion (2) requires the partial derivatives of the resid- 
uals with respect to the unknowns. These partial 
derivatives of the residuals are the coefficients in the 
tangential stiffness matrix K and are given by 


K a0 = d 2 II 
dq a dq@ 


/ 


C, 




dSj d£j 
dq a dqP 


dV + 



d 2 gj 

dq a dqP 


dV 

( 3 ) 


2 



The first and second integrals give the terms 
in the large displacement and geometric stiffness 
matrices, respectively. If the strains are equal to zero, 
the large displacement matrix obtained using the 
first integral and the nonlinear strain-displacement 
relations is identical to the matrix obtained by simply 
updating the coordinates and using the linear strain- 
displacement relations. If the strains are not equal to 
zero, there is a difference. This difference was verified 
numerically for 2-D elements. 

A Lagrangian formulation is used for NONLIN3D. 
For infinitesimal strain, the nonlinear strain- 
displacement relations (ref. 19) are 

£\ =u x + l/ 2 (u x u x + v x v x + w x w x ) 

£2 = Vy + 1/2 (UyUy + VyVy + WyWy) 

e?,-w z + 1/2 {u z u z + v z v z + w z w z ) . 

£4 = Uy + V X + ( U X Uy + V x Vy + W x Wy ) 

£5 —V Z +Wy + ( U Z Uy + V z Vy + W Z Wy) 

£6 = U Z + W X + {u z u x + v z v x + w z w x ) , 

where u, v, and w are displacements in the x-, 
y-, and ^-directions, and the subscripts x, y, and 
z indicate partial differentiation (e.g., u y = ^). 
The nodal values of u, v, and w are the unknowns 
referred to previously as q a . Note that £4, £5, and 
£6 are engineering shear strains. Since a Lagrangian 
formulation is used, the strains are based on the 
original configuration. For example, e\ is the axial 
strain of a line which was originally (i.e., before 
deformation) parallel to the re-axis. Although this 
line could be oriented parallel to the y-axis after 
deformation, the axial strain is still £1 and not £2. 

Substructuring 

The program NONLIN3D was designed to per- 
form analysis by substructures. A brief description 
of the substructuring technique is given here. More 
details can be found in reference 20. In addition to 
reducing computer memory requirements, substruc- 
turing allows the structure to be modeled as a combi- 
nation of linear and nonlinear components. For the 
configurations studied herein (fig. 1), linear analy- 
sis is appropriate everywhere except the majority of 
the postbuckled region. By substructuring into linear 
and nonlinear regions, expensive iterative solution is 
needed for only a fraction of the equations. 

In the current study, two substructures are used: 
one linear and one nonlinear. The procedure begins 
by obtaining a reduced stiffness matrix and load 
vector for the linear region. The reduced stiffness 
matrix can be treated as the stiffness matrix for 
just another type of element. Because of the large 


number of nodes, this element is referred to as a 
superelement. The stiffness matrix and load vector 
are “reduced” in the sense that only the nodes shared 
by the two substructures (the interface nodes) are 
included. 

The technique used for calculating the reduced 
stiffness coefficients utilized the formal definition of 
a stiffness coefficient: a stiffness coefficient is related 
to the restraint forces required to maintain unit 
displacement at one degree of freedom (DOF) and 
zero displacement at the remainder of the element 
DOF’s. Suppose there are to be n DOF’s in the 
superelement. These n DOF’s are restrained. One of 
these DOF’s is specified to have a unit displacement 
(and there are no other loads) and the governing 
equations are solved. The restraint forces at all n 
DOF’s constitute one column of the reduced stiffness 
matrix. This procedure is repeated for all n DOF’s. 
The reduced load vector is obtained in a similar 
fashion. All n DOF’s in the superelement are still 
restrained. However, the specified loads for the linear 
region are now applied. The reduced load vector is 
equal to the negative of the restraint forces at the 
restrained nodes. 

Once the superelement stiffness matrix and load 
vector are calculated, the analysis proceeds to the 
nonlinear substructure. Whenever the nonlinear 
stiffness matrix and load vector are formed, the in- 
teraction with the linear substructure is included by 
simply adding the superelement stiffness matrix and 
load vector. When the internally generated nodal 
forces are calculated to determine residuals, the con- 
tribution of the linear substructure consists of the 
product of the superelement stiffness matrix and the 
superelement nodal displacements. 

For the configuration analyzed, the delamination 
front is within the linear substructure. Hence, fur- 
ther work is required even after obtaining a converged 
solution for the nonlinear substructure. After ob- 
taining a converged solution, the displacements for 
the interface nodes are known. These displacements 
fully account for the effect of the nonlinear substruc- 
ture on the linear substructure. That is, the displace- 
ments in the linear substructure can be determined 
as though there were no other substructure, except 
that the magnitudes of the displacements at the in- 
terface nodes are specified. To reduce the computer 
resource requirements, it is usually advantageous to 
obtain multiple solutions for the nonlinear substruc- 
ture and then obtain multiple solutions for the linear 
substructure. 

Strain-Energy Release Rate Calculation 

The well-known virtual crack closure technique 
(ref. 21) served as the basis of the strain-energy 
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release rate calculation. This procedure determines 
Gj, Gfi, and Gju from the energy required to close 
the delamination over a short distance A a. The 
closure energy involves products of delamination 
front nodal forces and relative displacements be- 
hind the delamination front. The delamination front 
nodal forces can be determined by actually closing 
the delamination over Ao. Another technique, which 
requires only a single solution, assumes that the cur- 
rent delamination front nodal forces are the same as 
they would be if the delamination length were re- 
duced by Aa. The single-solution method was used 
herein. 

Figure 2 is a schematic of the delamination front 
region. The nodes of interest for the strain-energy 
release rate calculations are indicated by the filled 
circles. Because it is not appropriate to close the de- 
lamination over part of an element, there are four sets 
of nodes (indicated by the letters a, 6, c, and d ) which 
are used to calculate the closure energies. The rel- 
ative displacements are obtained by subtracting the 
displacements at nodes a' t and b\ from the displace- 
ments at nodes a % and b t , respectively. The forces are 
equal to the nodal forces transmitted across the de- 
lamination plane at nodes Cj and d % . These forces are 

obtained by evaluating the integral / Qy dV 
for all elements which are connected to nodes c t - or 
di and whose centroids lie above the delamination 
plane. There are two sets of energy products. One 
of the sets of energy products consists of the relative 
displacements for nodes a 2 - and a'- multiplied by the 
forces for nodes C{. The other set of energy products 
consists of the relative displacements for nodes b{ and 
b'l multiplied by the forces for nodes d{. The energies 
equal half of these products. 

Strain-energy release rate is a measure of energy 
per unit area. Hence, the energy products must be 
normalized by the appropriate areas. Unfortunately, 
there is not a simple exact way to determine the 
appropriate areas. The primary complication is that 
the midside nodes and corner nodes are “weighted” 
differently by the assumed element shape functions. 
The result is that, even if the strain-energy release 
rates were actually constant along the delamination 
front, there would be much larger energy products 
for the midside nodes than for the corner nodes. For 
example, in figure 2 the energy products associated 
with nodes c<i and c\ would be much larger than 
for those associated with nodes ci and C3. An 
approximate solution to this dilemma is as follows. 
The strain-energy release rate is not calculated for 
locations like C2 and C4 along the delamination front. 
Instead, the energy products associated with these 
locations are split evenly between the adjacent nodes. 


For example, the energy associated with location C3 
along the delamination front becomes 


^ ^0,30,303 "t" -^62^2^2 2 (■^'fl2 a 2 c 2 ^0484*24) 

(5) 

where E denotes the energy products associated with 
Gj, Gjj, and Gjjj, and the subscripts indicate 
the nodes involved. The area is approximated by 
the product of A a times the distance between the 
midside nodes on either side of the corner node being 
considered. For example, the area for node C3 is Aa 
times the distance from node C2 to node C4. 

If the delamination front is not parallel to one 
of the coordinate axes, it is preferable to add a 
coordinate transformation to the procedure outlined 
above. In particular, a local coordinate system is 
defined for each node along the physical delamination 
front (i.e., the Cj nodes). Figure 3 is a schematic 
of a delamination plane and the global ( xy ) and 
local ( x'y ') coordinate systems. This local coordinate 
system has one axis tangent to the delamination 
front, one axis normal to the delamination front, and 
one axis normal to the delamination plane. For all 
the cases considered, z and z' were parallel. The 
transformed nodal forces F x i , F y i , and F z i are defined 
as 


F x i = F x cos 6 + F y sin 0 
F y i = -F x sin 0 + F y cos 0 > 
F z > = F z 


( 6 ) 


The relative displacements are transformed simi- 
larly. The transformed forces and relative displace- 
ments are then used to calculate the energy prod- 
ucts. Figure 3 also defines the perimeter coordinate 
S, which is the distance along the delamination front 
measured from the y-axis. 

The procedure just outlined was implemented in 
two slightly different ways for the results presented 
here. The difference was in the way the nodal forces 
were calculated. Initially in the study, the nonlin- 
ear strain-displacement relations were used to calcu- 
late the nodal forces from the nodal displacements 
in the linear region. This is inconsistent, but if the 
region assumed to be linear is “exactly linear,” it 
makes no difference. The results for the mesh conver- 
gence study were obtained using this procedure. All 
the other results were obtained by using the linear 
strain-displacement relations to calculate the nodal 
forces from the nodal displacements in the linear re- 
gion. Because the linear region is not exactly linear, 
there is a difference in the results obtained using the 
two methods. The configuration used for the mesh 
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convergence study was also used in the parametric 
study, so results appear for that configuration using 
both methods. The second method is recommended 
by the author. 

Mesh Generation and Typical Model 

Figure 4 is an outline of the procedure used for 
mesh generation. This procedure is based on the 
procedure in reference 22. A two-dimensional model 
is swept through a 90° arc to generate a cylindrical 
3-D mesh. The outer part of the cylindrical mesh is 
then transformed to obtain a square boundary. Then 
an elliptical transformation is applied to obtain an 
elliptical delamination front. If the ellipse is longer 
in the y-direction than in the z-direction (i.e., b > a), 
the conformal transformation is 


x = x 


y' = y\ 1 + 


b 2 - a 2 
x 2 + y 2 


( 7 ) 


If a > b, the transformation is the same except 
that x and y are interchanged. To avoid a singularity 
in equation (7), nodes at zero radius were shifted to 
lie on an arc of very small radius, about 10~ 9 m. 

The transformation in equation (7) maintains the 
orthogonality of lines which were orthogonal in the 
modified cylindrical mesh (fig. 4) . This orthogonality 
at the delamination front simplifies the pairing of 
nodal forces and relative displacements in the strain- 
energy release-rate calculation. 

A peculiarity of the transformation in equa- 
tion (7) is the unusually close spacing of the elements 
close to the delamination front on the long axis of the 
ellipse. Figures 5(a) and 5(b) show the mesh before 
and after the transformation. Also, there appears to 
be a triangular element in figure 5(c). This element 
has two sides which are essentially collinear. A mod- 
ified transformation which results in more even mesh 
refinement is obtained by introducing a scale factor 
for the “stretching” as follows: 


y' = V 



b2 ~ “ 2 . F 
t 2 4- i/2 


F + 1 


where 


( 8 ) 


By making the parameter r a little less than the 
radius of the delamination front in the cylindrical 
mesh, orthogonality is maintained in the neighbor- 
hood of the delamination front during the elliptical 
transformation. 

After the elliptical transformation, the midside 
nodes are no longer at the middle of an element edge. 
Therefore, the coordinates of the midside nodes are 
recalculated as the average of the coordinates of the 
adjacent corner nodes. 

Figure 6 shows a typical finite element model. 
The elements are 20-node isoparametric hexahe- 
drons. Because of symmetry it is sufficient to model 
only one fourth of the specimen and impose the con- 
straints u = 0 on x = 0 and v = 0 on y = 0. There is 
also a constraint w — 0 on z = 0. This constraint was 
imposed to remove global bending from the analysis. 
In reality, there might be global bending (particu- 
larly if the buckled region is thick), but the amount 
of global bending depends on the region modeled and 
the boundary conditions at the external boundaries. 
Imposition of w = 0 on z — 0 simply removes overall 
specimen size and external boundary conditions as 
parameters to be considered in this study. The con- 
straint on w represents a laminate which is well con- 
strained globally. Of course, the imposition w = 0 
on z = 0 could also be viewed as an indication of 
symmetry about the z = 0 plane. This view implies 
the presence of two delaminations. 

Along the boundary x = W, all u displacements 
are specified to equal We x , where e x is the specified 
compressive axial strain. To initiate transverse de- 
flections, a transverse load was applied at the center 
of the delaminated region. After a converged solution 
was obtained, the load was removed, and solutions 
were obtained with only compression loading. 

Figure 7 shows a typical model after division into 
substructures. Most of the postbuckled region is in- 
cluded in the nonlinear substructure. The distance 
between the delamination front and the beginning of 
the nonlinear substructure was l. In all cases, l was 
approximately equal to the sublaminate thickness h. 
To check the validity of the substructuring, a crude 
model was analyzed with and without substructur- 
ing. The strain-energy release rates were negligibly 
different. For all the cases considered, the sublami- 
nate thickness h and base laminate thickness H were 
0.4 and 4 mm, respectively. 


F = ^ x2 + V 2 for r > Jx 2 + y 2 
r v 

and 

F = 1 for r < y x 2 + y 2 


Material Properties 

The objective of this paper is to consider only 
the effect of geometric parameters on Gy, G//, and 
Gm. Consequently, material properties were chosen 
to minimize the effect of material properties on the 
variations in G/, G//, and G ///. For quasi-isotropic 
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laminates, the in-plane stiffness is independent of the 
orientation in the xy plane. However, even for quasi- 
isotropic laminates, the flexural stiffness varies with 
orientation. Hence, even if the postbuckled region 
consisted of a quasi-isotropic group of plies, varia- 
tions would be expected in the strain-energy release 
rate along the delamination front solely because of 
the variation in flexural stiffness. Also, the proper- 
ties of the plies on either side of the delamination 
would be expected to at least affect the percentages 
of Gj, Gil-, and Gm . 

The simplified material properties chosen for this 
study are those for a “homogeneous quasi-isotropic 
laminate” throughout the entire specimen (buckled 
and unbuckled regions). These properties C are 
obtained as follows: 

5, = T(C S )‘ (9) 

k = 1 

where (Cy) fc are the constitutive properties for 
the ply in the 8-ply quasi-isotropic laminate 
(±45/0/90) s . With these properties throughout, 
there are obviously no stacking sequence effects and 
there is no variation of material properties with ori- 
entation. 

The lamina properties used to calculate Gy were 
selected to be typical of graphite/epoxy (ref. 23). 
The moduli and Poisson’s ratios are as follows: 

E n = 134 GPa 

E 22 = £33 = 10.2 GPa 

G \2 = Gn = 5.52 GPa 

G 23 = 3.43 GPa 1 j 

*T2 = ^13 = 0.3 

^23 = 0.49 


and approximate for large deflections. The equation 
for the central deflection w 0 is 

? — (?) 3 = o - 2 ^ ( 11 ) 

Finite element analyses were performed using a 
mesh similar to that in figure 6, but with a circular 
debond. The thickness and radius of the debonded 
region were 0.4 mm and 15 mm, respectively. The 
Young’s modulus E was 207 GPa, and the Poisson’s 
ratio was 0.3. 

Figure 8(a) compares the deflection at the center 
obtained from the closed-form solution with that of 
NONLIN3D. The agreement between the two anal- 
yses is excellent in the linear and initial nonlinear 
regions. There is a 10-percent difference in the two 
curves at the highest load level considered. This dif- 
ference is not surprising, since the closed-form anal- 
ysis is not exact for large deflections. 

The closed-form solution in reference 24 can be 
used to calculate the total strain-energy release rate. 
Because the configuration is axisymmetric, Gy is 
constant around the boundary. The strain-energy 
release rate is —(dU/da)/(2-iTa). The strain energy 
U can be calculated as the work done by the applied 
load (since all the work is stored as strain energy) as 
follows: 


U 


-I 

Equation (12) yields 


P dw 0 


U : 


2 . 30 ^+ 0 . 510 ^ 


(12) 


(13) 


a“ a * 

An expression for Gy is obtained by differentiat- 
ing equation (13) with respect to a and dividing by 
the circumference as follows: 


Results and Discussion 

The results of the analysis evaluation and para- 
metric study are described in this section. The analy- 
sis evaluation has several parts. First, the results for 
a transversely loaded laminate are discussed. This 
configuration was selected as a check case because 
it exhibits geometric nonlinearity and a closed-form 
plate solution is available. Then the results of the 
mesh convergence and allowable residual studies are 
discussed. Finally, the results of a limited parametric 
study are discussed. 

Transversely Loaded Plate 

A closed-form solution for a circular isotropic 
plate subjected to a central point load is given in ref- 
erence 24. This solution is exact for linear deflections 


Figure 8(b) shows Gy plotted against load. Re- 
sults are shown for the linear closed-form solu- 
tion (eq. (14) without the fourth-order w 0 term), 
the nonlinear closed-form solution (eq. (14)), and 
NONLIN3D. The results are plotted with log-log axes 
because of the wide range of the parameters. In the 
linear range, all three solutions agree very well. Even 
after nonlinear effects become important, the nonlin- 
ear closed-form solution and NONLIN3D agree very 
well. 

These comparisons indicate that NONLIN3D 
does account for geometric nonlinearity and that the 
technique for the strain-energy release-rate calcula- 
tion is valid. 


6 



Mesh Refinement 

Several configurations were analyzed as part of 
this study. Consequently, it was neither practical 
nor warranted to perform a convergence study for 
all cases. Instead, a systematic convergence study 
was performed for a single configuration. Mesh 
refinements for other configurations were selected 
based on the results of the convergence study. 

The configuration selected for the convergence 
study had a circular delamination with a radius of 
15 mm. The sublaminate thickness h was 0.4 mm 
and the base laminate thickness H was 4 mm. The 
overall laminate width W was 50 mm. The ma- 
terial properties were those for the homogeneous 
quasi-isotropic laminate described previously. Fig- 
ure 9 shows the extremes of refinement used for the 
2-D meshes. The elements in the coarse mesh were 
subdivided to obtain the refined mesh. For the coarse 
mesh, only two elements were used to model most 
of the buckled region. As described previously, the 
2-D meshes were swept through a 90° arc to generate 
the 3-D meshes. Figure 10 shows four of the meshes 
generated. As shown in the figure, the number of 
slices of elements was also varied. Models with 4 and 
8 slices are shown; a 12-slice model was also used. In- 
formation on all the models used in the convergence 
study is given in table 1. As shown by the figures 
and the table, a fairly wide range of refinement was 
examined. 

Since strain-energy release rates were of primary 
importance, variations in Gj, Gjj, and 6’/// were 
used to determine the adequacy of the mesh refine- 
ment. Figure 11 shows the distribution of Gj and 
Gj i along the delamination front for three strain lev- 
els and four models (models 1, 3, 4, and 6). Only 
symbols are shown for the crudest mesh, model 3. 
These meshes bracket the entire range of refinement 
in table 1. The mode III component Gjjj was negli- 
gible for all cases. Of interest here are the differences 
in the results obtained using the various meshes. Ex- 
cept for model 3, which only had 4 slices of elements, 
the results from all the models are essentially equal. 
Even a coarse 4-slice model gives the correct trends. 
Apparently, a fairly crude model is sufficient to cal- 
culate Gj and Gjj. The significance of the magni- 
tude and the distribution of Gj and Gjj is discussed 
subsequently. 

Allowable Residual 

The program NONLIN3D solves the governing 
nonlinear equations by using a modified Newton- 
Raphson solution procedure. During each iteration, 
the internally generated forces are compared with 
the externally applied loads. The differences are the 
residuals. If all the residuals are identically zero, the 


governing equations are exactly satisfied. Of course, 
in practice, exact agreement is seldom obtained. It- 
eration could continue until the algorithm’s best ap- 
proximation of zero is obtained. The size of this “nu- 
merical” zero depends on the computer and the type 
specifications in the program. However, negligible 
residuals are, in general, orders of magnitude larger 
than a numerical zero. 

To determine what is a negligible residual, three 
residual tolerances were considered: 1000 N, 1 N, 
and 0.0001 N. A laminate with a 30 x 60 mm 
delamination was analyzed for five strain levels. The 
range of strains was such that the maximum lateral 
deflection for the buckled region varied from about 
0.6 to 2.2 times the thickness of the buckled region. 
The tolerance of 1000 N gave erroneous results. The 
other two tolerances gave virtually identical results 
except for the lowest strain level ( e x = —0.001), 
for which there were differences of about 6 percent. 
Figure 12 shows Gj and Gjj for e x = —0.001 (i.e., 
the worst case). Results for the other strains are not 
shown, since the differences are very small. 

Based on these results, a tolerance of 0.0001 N 
was selected for all the analyses. A somewhat larger 
residual could probably have been tolerated. How- 
ever, the residuals tend to decrease quite rapidly dur- 
ing the iterations, so there is little to be gained (in 
terms of reduced cost) by trying to specify the largest 
acceptable residual. Also, a larger tolerance might 
give poor results for some cases in which the loads 
are small. 

Parametric Study 

The parameters considered were strain level, de- 
lamination shape, and delamination size. First, the 
deformation of the postbuckled region is discussed. 
Then the calculated strain-energy release rates are 
presented. 

Figure 13 shows lateral displacement in the mid- 
dle of the delaminated region plotted against ax- 
ial strain for two circular and two elliptical delam- 
inations. The dimensions of the delaminations are 
shown in the figure. The curves indicate that when 
the buckling load is exceeded, the displacement in- 
creases rapidly at first with increased strain. Then 
the rate of increase in displacement decreases. Obvi- 
ously, the response is quite nonlinear. 

Figure 14 shows deformed finite element meshes 
for a circular and an elliptical delamination. The 
displacements have been multiplied by 10 to improve 
visualization. The deformed shape is relatively sim- 
ple except near the delamination front. For both 
cases, the delamination front is open near the y-axis. 
However, for the circular delamination, the delam- 
ination faces actually overlap near the a;-axis (i.e., 
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near the y = 0 plane). Even at strains smaller and 
larger than those shown, the circular delamination 
exhibited overlapping. For the elliptical delamina- 
tion, “small” strains result in opening of the delami- 
nation along the entire length. At larger strains (not 
shown), the elliptical delamination also exhibited 
overlapping. Strictly speaking, constraints should 
be added to prevent overlapping of the delamina- 
tion faces. However, including constraints to prevent 
overlapping further complicates an already compli- 
cated stress analysis problem. Consequently, no con- 
tact constraints were added for any of the results 
presented in this paper. In the results that follow 
(figs. 15 to 17), dashed lines are used for the Gj and 
Gjj distribution curves in regions where overlap oc- 
curred. 

Figure 15 shows the Gj and Gjj distributions for 
a circular delamination for five strain levels. This 
is the same configuration used for the convergence 
study. The G/// distribution was essentially zero for 
this and all other cases considered in this study. In 
general, Gm would not necessarily be expected to 
be zero. The strain-energy release rates are plotted 
in figure 15 using the perimeter coordinate 5. This 
coordinate is zero where the delamination front meets 
the y-axis and is maximum where the delamination 
front meets the z-axis. Both G/ and G// show 
large variations along the front and are largest at 
5 = 0. There is overlapping of the delamination 
surfaces over a large portion of the front, as indicated 
in figure 15(a). Although G/ is larger than G// 
for the five strain levels, the difference is not large; 
this is definitely a mixed-mode situation. Since both 
G i and G// are largest at S ~ 0 , delamination 
growth would be expected to occur preferentially 
perpendicular to the load direction, that is, in the 
^-direction. 

Since a circular delamination is expected to be- 
come elongated perpendicular to the load direction, 
a 30 x 60 mm elliptical delamination was ana- 
lyzed. Figure 16 shows the distribution of G/, G//, 
and Gt for this elliptical delamination. There is 
a large variation of both Gj and Gjj along the 
front. Note that the location of the G/ peak shifts 
slightly with strain level. In contrast to the circu- 
lar delamination, the peak values of G/ and Gjj 
occur at different locations. Also, the peak value 
of G// is larger than the peak value of G/, ex- 
cept for the case e x = —0.005. The total strain- 
energy release rate (fig. 16(c)) varies significantly 
along the front for the larger strains, but for the 
smaller strains the magnitude is almost constant. 
Depending on the choice of growth criterion, very 
different predictions of even the direction of growth 
are possible. A criterion based only on G/ would 


predict growth perpendicular, or nearly perpendicu- 
lar, to the load direction for all the strain levels. A 
criterion based only on G/j would predict growth 
parallel to the load direction. For the smallest 
strains, a criterion based on total strain-energy re- 
lease rate would predict almost uniform growth along 
the delamination front. 

Figure 17 shows G/ and G// for a 60 x 60 mm 
delamination. Comparison of figures 15 and 17 shows 
the effect of delamination size on Gj and G// for 
two circular delaminations. Both delaminations were 
subjected to the same strain levels. Figures 15(a) 
and 17(a) show that the larger delamination is closed 
(actually overlapping) over more of the delamination 
front. Also, the distribution in the overlapping region 
is more complicated for the larger delamination. This 
is because the strains for the larger delamination are 
larger multiples of the bifurcation buckling strain. 
This conclusion was verified by subjecting the smaller 
delamination to higher strains. (These results are not 
presented in this report.) The larger delamination 
has a much larger G/ for the region near 5 = 0. 
Figures 15(b) and 17(b) show that G// is also larger 
for the larger delamination near 5 = 0. Hence, 
unstable extension of the delamination might be 
expected once it begins to grow. However, based on 
the calculated strain-energy release rates, a circular 
delamination is not expected to grow self-similarly 
into a larger circular delamination. It should become 
elliptical. 

Figures 15 and 16 show results for 30 x 30 and 
30 x 60 mm delaminations. Comparison of figures 15 
and 16 shows that, except for the lowest strains, the 
smaller (circular) delamination actually has a larger 
G/. Hence, the growth rate based on G/ is expected 
to be larger for the smaller delamination. As pointed 
out previously, the distributions of Gjj for circular 
and elliptical delaminations are quite different. Even 
the location of maximum G// is different. The peak 
values of G// for the two delaminations are similar. 
Therefore, based on Gjj, the growth rate should 
be about the same for both delaminations, but the 
direction of growth would be different. 

Conclusions 

A three-dimensional, geometrically nonlinear, 
finite element program, NONLIN3D, was developed 
to calculate the strain-energy release rate compo- 
nents for a postbuckled embedded delamination. The 
program NONLIN3D uses substructuring to reduce 
the number of equations which must be solved iter- 
atively. A new technique was developed for calcu- 
lating strain-energy release rates G/, Gjj, and G///. 
This technique is based on virtual crack closure. Var- 
ious checks of NONLIN3D and the finite element 
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models were performed to establish confidence in the 
results presented herein. A parametric study was 
performed for postbuckled elliptical delaminations 
in a “homogeneous quasi-isotropic laminate.” Based 
on the parametric analysis, the following conclusions 
were reached: 

1. The problem is definitely mixed-mode. In 
some cases, Gj is larger than Gjf, for other 
cases, the opposite is true. 

2. In general, there is a large gradient in the 
strain-energy release rates along the delami- 
nation front. 

3. The locations of maximum Gj and G;j de- 
pend on the delamination shape and the 
applied strain. 

4. The mode III component was negligible for all 
cases considered. 

5. The analysis predicted that parts of the 
delamination would overlap. The analysis pre- 
sented herein does not impose contact con- 
straints to prevent this overlapping. Fur- 
ther work is needed to determine the effect of 
allowing the overlapping. 

NASA Langley Research Center 
Hampton, Virginia 23665-5225 
May 26, 1988 
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Table 1. Statistics on Models Used in Convergence Study 
[Substructures 1 and 2 are linear and nonlinear, respectively] 


Model 

Number of nodes/ 

Slumber of elements 

Number 

2-D model 

3-D model 

Substructure 1 

Substructure 2 

1 

143/36 

1719/288 

1475/256 

287/32 

8 

2 

153/38 

1841/304 

1475/256 

409/48 

8 

3 

173/42 

1129/168 

799/128 

353/40 

4 

4 

173/42 

2085/336 

1475/256 

653/80 

8 

5 

173/42 

3041/504 

2151/384 

953/120 

12 

6 

533/152 

6325/1216 

5173/1024 

1221/192 

8 
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s= o 

9 = 0 


S = Length of 
delamination front 

Figure 3. Global and local coordinate systems and perimeter coordinate S. 



2-D mesh 




Elliptic mesh 



Modified cylindrical mesh 


Figure 4. Procedure of generating three-dimensional finite element models. 



(a) Before transformation. 


(b) After elliptic transformation. 



(c) Close-up of region indicated 
by arrow in (b). 


(d) After modified elliptic 
transformation. 


Figure 5. Transformed meshes using elliptic transformation of reference 22 and the modified transformation. 
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Figure 6. Typical finite element model. 

Nonlinear substructure Linear substructure 



Delamination front — ' ' \ \ ' 

Figure 7. Typical division of model into substructures. 
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(b) Total strain-energy release rate. 

Figure 8. Analysis of transversely loaded plate with circular delamination. 



(b) Refined 2-D model which has 533 nodes. 

Figure 9. Range of two-dimensional mesh refinement for convergence study. 



(c) Model 4. (d) Model 6. 

Figure 10. Several of the three-dimensional meshes used in convergence study. 
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(b) Mode II. 

Figure 11. Strain-energy release rates calculated using different meshes. 
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(a) Mode I. 
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(b) Mode II. 

Figure 15. Strain-energy release rate distributions for a circular delamination. (Diameter = 30 mm.) Dashed 
lines indicate overlapping. 
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(b) Mode II. 

Figure 16. Strain-energy release rate distributions for an elliptical delamination. (2a x 26 = 30 x 60 mm .) 
Dashed lines indicate overlapping. 
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(b) Mode II. 

Figure 17. Strain-energy release rate distributions for a circular delamination. (Diameter = 60 mm.) Dashed 
lines indicate overlapping. 
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